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Abstract 



Oh 

< 

We introduce an iterative method for computing the first eigenpair (A p , e p ) for the p- 

SLaplacian operator with homogeneous Dirichlet data as the limit of (fj, q u q ) as q — > p~, where 
i— —i u q is the positive solution of the sublinear Lane-Emden equation —A p u q = fi q u q with the 

same boundary data. The method is shown to work for any smooth, bounded domain. 
Solutions to the Lane-Emden problem are obtained through inverse iteration of a super- 
solution which is derived from the solution to the torsional creep problem. Convergence of 
u q to e p is in the C 1 -norm and the rate of convergence of \i q to X p is at least O (p — q). 



£C) Numerical evidence is presented. 

t-h Keywords: p-Laplacian, first eigenvalue and eigenfunction, inverse iteration, Lane-Emden problem, 
^ torsional creep problem. 



X 1 Introduction 



In this paper we develop an iterative method to obtain the first eigenpair (A p , e p ) of the eigenvalue 
problem 

-A p u = A \u\ p ~ 2 u in Q, 



u = on <9f2, 



(1) 



*E-mail addresses: rodney@mat.ufmg.br (R. J. Biezuner), jed@59A2.org (J. Brown), grey@mat.ufmg.br (G. 
Ercole), eder@iceb.ufop.br (E. Martins). 
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where A p u := div (|Vm| p ~ 2 Vtt), p > 1, is the p-Laplacian operator and C M. N , N ^ 2, is any 
smooth, bounded domain. The p-Laplacian equation appears in several mathematical models 
in fluid dynamics, such as in the modelling of non-Newtonian fluids and glaciology [SI EH [2IH 
HI], turbulent flows [22], climatology [21] nonlinear diffusion (where it is called the iV-diffusion 
equation; see [12] for the original article and [28] for some current developments), flow through 
porous media [IS], power law materials [6] and in the study of torsional creep [31] . 
The first eigenvalue X p of Q is variationally characterized by 

X p = min R (u) > 

uew£'*(n)/{o} 

where i? is the Rayleigh quotient 



in 


Vu 


\ p dx 


hi 


\u[ 


p dx 



The first eigenfunction e p of <Q is characterized by the fact that the minimum of R is attained 
at e p , so that 

A = / n |Ve P r<fa 

It is well-known that X p is isolated and simple, and that the corresponding eigenfunction e p G 
C 1,a (Q) can be taken positive. Since R is homogeneous, we may assume HepHoo = 1, where || • 
stands for the L°°-norm. 

In the one- dimensional case the first eigenpair (X p , e p ) is explicitly determined by solving 
the corresponding ODE boundary value problem. If Q = (a, b), then X p = (ix p / (b — a)) p 1 and 
e p = (p — l) _1 ^ p sin p (7T p (x — a) / (b — a)), where tc p := 2(p — I) 1 /* 5 J Q (1 — s p )~ l / p ds and sin p is a 
27T p -periodic function that generalizes the classical sine function (see [TTj HO])- 

When p = 2, we have A p = A, the Laplacian operator, whose first eigenpair (X p , e p ) is well- 
known for domains with simple geometry (that is, domains which admit some kind of symmetry); 
for more general domains it can be determined by several numerical methods (see [121 EHl EH 
[30l [36]). However, if p ^ 2 and N ^ 2, the first eigenpair is not explicitly known even for simple 
symmetric domains such as a square or a ball, and there are few available numerical methods to 
deal directly with the eigenproblem ([!]) in these domains (see (TDJ [EH HSJ dSl EBJ HZ])- 

On the other hand, several numerical methods are available to solve homogeneous Dirichlet 
problems for the (Poisson) p-Laplacian equation in the form 

—A p u = f (x) 

when / depends only on x G Q (see [21 El El 1231 12SJ 115] )• This fact motivated the development 
of an inverse iteration method by some of the authors for finding the first eigenpair in [TU]. If Q 
is a N- dimensional ball, the convergence of the method was established and numerical evidence 
for its applicability when Q is a 2-dimensional square were also presented. In the special case 
of the Laplacian operator, the method was proved to work in general domains and can be also 
used to obtain other eigenpairs (see [12])- However, since the method was based on the iteration 
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of the nonlinear p-Laplacian equation in Q, the difficulties in dealing with the nonlinearity on 
the right-hand side of the equation prevented showing that the method works in any domain and 
any p > 1 . 

In this work we consider a different inverse iteration approach, also based on the solution 
of the Poisson p-Laplacian equation, but built around an eigenproblem which has a sublinear 
nonlinearity on its right-hand side. This type on nonlinearity is more manageable and we are 
able to prove that the iterative method works for any smooth, bounded domain. It is based on 
obtaining positive solutions v M for the Lane-Emden type problem 

-A p v = /i \v\ q ~ 2 v in Q, ,„s 

= on on. u 

After rescaling, /i and v^ q produce a family of pairs {(/i g , u q)}i <q<p converging to the first eigen- 
pair (A p , e p ) when q — > p~, the convergence u q — > e p being in C 1 (Q). We will now describe the 
method in more detail. 

It is well known that for each fixed /i > 0, problem ^ has a unique solution v^g, if 1 < q < p 
(see [32J ) . If q = p, we have the p-Laplacian eigenvalue problem. If q > p, positive solutions of 
(J5]) usually are not unique. A nonuniqueness result for ring-shaped domains is given in [7j when 
q is close to the Sobolev critical exponent p* (p* = Np/ (N — p), if 1 < p < N, and p* = oo, 
if p ^ N). On the other hand, as proved in pQ, positive solutions are unique when Q is a ball, 
while for general bounded domains the uniqueness of positive solutions that reach the minimum 
energy (ground states) was established in [24J under the conditions 1 < p < N and I < q < p*. 

Now, in order to construct the approximating sequence to the first eigenpair, first choose any 
/i > and a sequence (g n ), 1 < q n < p, such that q n — > p~ . It is important to notice that /i 
need not to be taken close to X p . This point is crucial, since good a priori estimates for X p are 
hard to obtain. For each q n we need to solve the Lane-Emden problem (|2]) in order to find v Mn , 
which is a degenerate nonlinear problem almost as hard to solve as the eigenvalue problem for 
the p-Laplacian ([!]) itself. In order to obtain the solutions v^ qn we first solve the much easier 
torsional creep problem 

-A p =1 in Q, , . 

i> = on dQ. 1 } 



Then compute k p = ||0||^ o p and set 



(4) 



0o is a supersolution to (J2|. One immediately sees that the easiest choice is /i = k p , so that 
0o — 0/||0||oo- Now apply an inverse iteration to O > finding a sequence of iterates ((f> m ) which 
satisfy 

-A p( f) m+1 = ^(j)^ 1 in VI, 
Pm+i = on dQ. 

This can be done by a number of numerical methods. Finite volume based methods are presented 
in [H [25] ; finite element based methods are also available (see [31] and the references therein). 
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After a pre-established tolerance limit has been reached at some <p m , where m is a function of /i 
and q n , set 



and define u qn and // 9n as 



// 



^ aIld U «» 



y M>9n II OO 



^iQn II oo 



In Theorem 
value for q c 



8 we show that jj, qn — > \ p and tt 3n — > e p in C 1 (fi) when q n — » p . Choosing a 
ose to p will give an approximation for the first eigenpair of the p-Laplacian. The 



procedure is summarized in Algorithm 1 below. 



Algorithm 1 Inverse iteration for the first p-Laplacian eigenpair (X p , e p ) 



set [i 
set q 

solve — A r 



1 in f2. 



on dQ 



set 4>q = (fi/k p )p-i dp 
for m = 0, 1, 2, . . . do 

solve -A p m+ i = /i(/>m 1 in 0, 
end for 

return fj,/ ||0 m +i||^" 9 
return </> m+ i/ ||0 m +l IL 



(an arbitrary positive number) 
{q should be chosen close to p) 
(torsion function) 

(supersolution) 

f) m+ i = on (90 (Inverse iteration sequence) 

(first eigenvalue X p ) 
(first eigenfunction e p ) 



However, we are able to produce a much more robust algorithm which is also easier to apply in 
practice. In Algorithm 2 below one does not need to use an arbitrary parameter /i nor to compute 
the value of the constant k p . Normalization of the iterates at each step increases robustness and 
thus it should be the algorithm of choice. 



Algorithm 2 Inverse iteration for the first p-Laplacian eigenpair (A p , e p ) with normalization 

1 

2 
3 
4 
5 
G 
7 



set q 

solve —A p (pQ = 1 in 0, < 
for m = 0, 1, 2, . . . do 

solve -A p (f) m+ i = (<p m 
end for 

return 1/ H^+ill^ 1 
return <f> m+ i/ H^WiH^ 



/ 



= on <9Q 

^mlloo) 9 " 1 in n, 



(q should be chosen close to p) 
(torsion function) 

on d£l (Inverse iteration sequence) 

(first eigenvalue X p ) 
(first eigenfunction e p ) 



The outline of the paper is as follows. In Section 2 we present some preliminary results 
that will be used in the sequel. The sequences of approximates for both algorithms are built in 
Section 3 and the proof of their convergence to the first eigenpair is given in Section 4. In Section 
[5] we present some numerical results for the unit ball of dimensions N = 2,3,4 using the first 
algorithm, and for the two-dimensional square and the three-dimensional cube and torus using 
the second algorithm. 
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The main advantage of the method presented here, besides its applicability to general domains, 
is that approximations to both \ p and e p are obtained with the desired precision by an iterative 
process which is numerically simple and, in the case of a ball, also explicit. 



2 Preliminary results 

In this section we state simple versions of some results on the p-Laplacian. We begin with the 
following comparison principle (see [TH] for a more general version). 

Lemma 1 For i G {1,2}, let hi G C (Jfj and U{ G W 1)P (VL) be such that —A p Ui = hi in Q. If 
h\ < h2 in Q and u± ^ U2 on dQ, then u\ ^ U2 in Q. 

The following result is a simple version of a general result proved in the classical paper [3H] 
of Lieberman. 

Theorem 2 f3R Thm 1] Suppose that u G W 1,p (Q) is a weak solution of the Dirichlet problem 

—A p u = f(x, u) in Q, 
u = on dQ 

where f is a continuous function such that 

1/(2,01 < A forall{x,0enx[-M,M) 

for positive constants A and M. 

V \\ u \\oo ^ then there exists < a < 1, depending only on A, p and the dimension N, 
such that u G C 1,Q (JTj ; moreover we have 

IMIci.«(n) ^ c> 

where C is a positive constant that depends only on A, p, N and M. 

Thus, denoting by is the solution of the torsional creep problem ^ in the domain f2, one 



can easily verify using (12) and the comparison principle in balls that < ^ M in Q for some 
positive constant M. Hence, Theorem 2 implies that G C 1,a (Cl) for some < a < 1. 
For the next lemma set 

kp ■= U\t P > 0. (6) 



Lemma 3 k p ^ X p . 



Proof. Let e p be the first eigenfuncion associated with X p satisfying ||e p || = 1 in Q. Since 



- A p e p = \peP~ 1 ^X p = -A p ^A^ 1 0J in fi, 
e p = = Ap _1 on dil, 
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it follows from the comparison principle that 



Hence, 



from what follows our claim. 



O<e p ^A^ 1 in a 



1 — IMloo ^ 1 II 1 'I! x ' 



Remark 4 It follows from Picone's identity (see that, in fact, the inequality is strict, that 
is, k p < X p (for details, see Lemma 8.1]). 

The following result is well-known and follows from Theorem [2] 

Theorem 5 Let — A p 1 : C l (fi) — > W Q ' P (O) be the operator defined as follows: for each v G 
C 1 (Q) let — 1 v :=m6 Wq' p (O) be the unique solution of the Dirichlet problem 

—A p u = v in £7, 
u = on dQ. 

Then —A' 1 is continuous and compact. Moreover, —A^v G C 1,a (fl) for each v G C 1 (Q) . 

In the remainder of the paper (A p , e p ) denotes the first eigenpair of ([!]), denotes the torsion 
function of Q and k p := . 



3 Construction of the sequence of approximates 

As mentioned before, if 1 < q < p, then for each fi > the Lane-Emden problem 

— A p v = fi \v \ 9 ~~ 2 v in a 
v = on <9a 



(7) 



has a unique positive solution v^g, which can be obtained via standard variational, and therefore 
non- constructive, arguments. The existence and uniqueness of solutions of ([7]) in the case 1 < 
q < p implies that the map fi i— > Vp >q is well-defined and monotone, in the sense that /i X < /i 2 

implies v^ uq < v^ q in VL, since v^ uq = (/ii//i 2 ) 1/(p_9) for an y A*i> /4a > 0. 
The basis of our constructive method is given by 

Theorem 6 Suppose 1 < q < p. For each fi > the unique positive solution v^ q G C 1 '" 

rl, P 



Wq' p (fl) of <m satisfies 



p / \ n p 



Moreover, v^ q is the limit, in the C 1 (Jfj norm, of the sequence 
iteratively defined by 



i 



k p 

and, for n ^ 1, 

-A p v n+1 = /jlvI* 1 in tt, 
v n+ \ = on dVt. 

Proof. DeBne := ,ne p and ,„ := where 



y i : = hr ra- ( 9 ) 



(10) 



m: =ur and M: Hc 

We have 

- A P ^,g < and - A p ?7 Mi5 ^ //v*" 1 in Q. (11) 

Indeed, in f2 in we have 

-A v = A ?; p_1 = A v p ~ q v q ~ 1 = A fmp "l^ 9 ?; 9-1 < A m v ~ q v q ~ x = nv q ~ l 

and ^ 

-ApU^ = k p M?- 1 = k v MV~ q M q ~ x > k p MP~ q (7^-) = fivft. 



Since v q = = v^ q on <9f2 the inequalities in (11) mean that v and v^ q are, respectively, sub- 
and super solutions for 0. 

Moreover, and are ordered, that is v„ q ^ v^q in fi. For, since k p ^ A p , we have 

\ p m p ~ l = \ p Upj P ~" 

q — 1 g — 1 p — 1 



p-1 / 1 \ P-9 p-1 / 1 \ P-Q ( U \ P-9 . 

= h£' =m * 



whence 



-A^ = V£? < \ p mv- 1 ^ k p M?- 1 = -A p v m 



in Q. Thus, since u„ = = on <9fi, we obtain u ^ in Q by applying the comparison 
principle. 

Since u h-> {iu q ~ x is increasing and y_ /iq ^ v^g in f2, the comparison principle also implies 
that the sequence {t> n } defined by the iterative process (10) starting with the supersolution v m 
satisfies 

v m < v n+1 ^v n ^ v m in fi. 
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Hence ■ y n converges to a function v m a.e. in Q. Since HtVilloo ^ ll^glloo = ^ f°ll° ws from 
Theorem m that {v n } C C l > a (H) for some < a < 1 (which does not depend on n) and that 



for some positive constant C which is independent of n. 

Thus, from Arzela-Ascoli theorem we conclude that v n — > v in the C 1 norm. 
Now, the continuity of the operator —A" 1 : C 1 (Q) — > Wq' p (Q) permits passing to the limit 



_ p 

in (10), which yields that v m G C 1 (fi) fl W ' p (fi) is a solution of (7| satisfying 



< v^ q ^ v^ q ^ v^ q in VI, 

proving (TsJ) . The regularity v^q G C 1,a (Jl) follows from Theorem [2J ■ 

This iterative process is also known as inverse iteration since v n+ i = — A p 1 (fiv^ 1 ). It is 
essentially the sub- and supersolution method starting with the supersolution v^ q ; the solution 
v^ q that it produces is characterized as the maximal solution between v and v^ q . 

If one starts the iteration with the subsolution then one obtains an increasing sequence con- 
verging to the minimal solution between v and Because of the uniqueness this minimal 
solution coincides with v^ q . However, in order to compute the minimal solution from this itera- 
tive process, it is necessary to know a priori a subsolution, which is exactly one of the unknowns 
that we wish to find by applying the method. 

On the other hand the supersolution v^q is easily obtainable since it involves the solution of 
the simpler problem 

For example, if Q = Br(xq), the ball centered at Xq G IR with radius R > 0, then it is easy 
to verify (see also below) that the torsion function is the radial function 

p — 1 / p p \ 

4> ( r ) = — ( R^ 1 — \r\ J , r = \x — Xq\ ^ R. (12) 



pN 



We then obtain that 



and 



K = U\t p = ^-J-^Y 1 (13) 



Rp \p — l 




where r = \x — Xq\ . 

In this case the sequence v n converging to v^ q is given recursively by the formula 

v n+ i (r) = f ( [ (^Y' 1 fjLvM^ds^) " 1 d9 (14) 



where v (r) = v m (r 



S 



This integral formula follows from the more general fact: the Poisson problem 

-A p u = f(\x\) in B R (x ) 
u = on OBr^xq) 

is equivalent to the ODE boundary value problem 

- (r N ^ \u'\ p - 2 u')' = r w "7(r), < r < R 
u'(0) = = u{R) 

for radial solutions u = u(r), r = \x — xq\ . Hence, after two integrations of the ODE taking into 
account the boundary conditions one obtains the following integral expression 

u(r) = J* (jf (J)" _1 /(>)<k) P ^ d0 ( 15 ) 

for the solution u = u(\x — xo\) of the Poisson problem. In particular, when f(r) = 1 this integral 
form might be simplified in order to find the expression (12) for the torsion function of Br. 

In our method, in order to compute the first eigenpair (X p ,e p ), we fix a positive value /x > 
and choose q close to p~ . Then, we apply the inverse iteration of Theorem [6] starting with the 
supersolution 



i 

p-q 



to obtain approximations for the function v^ q . Hence, 

t 1 -A V , 



— » X p and j. — ^j. > e p (in the C 1 norm) 



11 r H 7) 
\ v »,q\\oa II ^A*,? I loo 



as q — > p, a result that we prove in the next section. 

For the construction of the normalized sequence of Algorithm 2 one needs the following result: 

Theorem 7 Suppose 1 < q < p. Then the normalized sequence {w n / Halloo} where w n is defined 
by 



w := 1 and 



-A p w n+ i = { ,, Wn , — | in Vl 



Fnlloo 



w n+ \ = on dVt 

converges in the C l (Vt) norm to v q / WvgW^ where v q G C 1,a (Q) fl Wq ,p (Q) is the solution of |?|) 
with fi = k p . 



Proof. Let {v n } be the sequence defined by 

Vl - and 1 v „ +1 = o on an. ' ( ) 
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It follows from Theorem the sequence {v n } is decreasing and converges in the C 1 (r2) 

to the solution v q G C 1,a (Q) fl W Q ' P (Q) of the Lane-Emden problem 

— A p v = k v v q ~ x in Q 
v = on dQ 

Since u>i = we have that 

i_ 

w 2 = k p p ~ 1 v 2 

and, in particular - — \ — = - — \ — . In fact, this follows from the comparison principle: 

iFalL INL 



-A p (k p p - 1 v 2 ) = kZ 1 (-A p v 2 



K X K V \ 1 = ( TTTiT ) = ( I, W \ ) = -A p w 2 . 



Repeating this procedure we obtain 



-a p ( k p ihlU^ v 3 ) = K 1 Mi^-v (-a p v 3 

- K r> \\ v 2 no k p v : 



v 2 Ly v h L/ 



that is 



and 



1 q-l 

p— 1 II || p— 1 



W 3 = k p P 1 ||w 2 ||oo P 1 ^3 
W 3 W 3 



IKIL INL' 
Therefore, by an induction argument we conclude that 

i i-i 



P-i II || p-i .. . . i V n+1 



w n+1 = k p p \\v n \\ 00 p v n+ i and 



K+iH^ ||f n+ i| 

for all n ^ 2. Hence, it follows from Theorem [6] that 



Ihn+llL Ikn+llL ' Halloo 
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4 Convergence of the method 

Theorem 8 For \x > and for each 1 < q < p set 



u q := ] r^j r , (18) 



where v^ q G C l,a (fi) is the unique positive solution of |?]) ; and 



to ■= |i7-ipf • ^ 
II ^.alloo 



Then ji q — > \ p and u q — > e p in C 1 {ft) as q —> p . 

Proof. Since llttJI = 1 and 

ii y ii oo 



P 1 n np-l M,9 II ||P-? 1 t Jj 1 a q ' 



F/mIIoo II^V>9lloo 



we have that u q is the unique solution of the problem 



As a consequence of dsT) we have 



-A p u q = fi q u q q 1 in Q, 
u q = on dfl. 



— < \\v \\ p ~ q < — 



(20) 



< ( "A" e T < u, < f^" * i„Sl (21) 



and 

fc P < ^ < Ap. (22) 

Since 

< 1 < A p , 

it follows from Theorem [2] the existence of constants < a < 1 and C > independent of q such 
that u q e C 1,Q (H) and 

IWIcWli) ^ ^ for all 1 < g < p. 

Using the compactness of the immersion C 1,Q (Qj ^ C 1 (fi), letting g n — >■ p we get, up to a 
subsequence, /i 9n — >■ A 6 [fcp, A p ] and u gn — > u in C 1 Taking the limit in (20), we conclude 
from Theorem [5] that u must satisfy 

— A p u = Aw p_1 in Q, 
u = on <9fi, 

and IHloo = 1, whence A = X p and u = e p because A is an eigenvalue and u ^ is a corresponding 
eigenfuntion that does not change the signal in Q (note from (21) that u > in Q). Since these 
limits are always the same, that is, do not depend on particular subsequences, this ends the 
proof. ■ 
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Corollary 9 Using the notation of Theorem^ it follows that 



in the C 1 (r2) norm and 

Win I , 

loo ^ /v p 

as n —7- oo and q — > p~ . 



Proof. For each 1 < g < p let t> g denote the solution of the Lane-Emden problem (17). It follows 
from Theorem [8] that v q / H^H^ converges in the C 1 -norm to the first eigenfunction e p and that 



Using the notation of Theorem [7] we have 



lim lim Wn — = lim lim - — j — = lim - — j — = e p . 

g->f>- n->oo ky n q-*p- n-^oo \\ V n I g->p" % L 



Moreover, 

lim Hwn+ill^T 1 = lim k' 1 Ikn+ill™ 1 ||f n || 1-9 = K 1 IKIIfT 9 



and hence 



lim lim ||w7 n+ i||^ p = X p . 

q->p~ n^oo 



Next we prove an error estimate in the approximation of X p by p, q or, alternatively, by the 
scaled quotient 



\\v 

A q := H 



p 



where ||-|| r denotes the norm of the U (Q) , that is, ||w|| r = (j n \w\ r dx) r . 

The upper bound A q together with the lower bound \i q allows one to better control the 
accuracy of the approximation to X p . 

Theorem 10 There holds: 

(i) X p < A,. 

(ii) A q — » X p as q —¥ p~. 

(Hi) There exists a positive constant K which does not depend on q such that 

*C max {{X p - fM q ) , (A, - X p )} < K (p - q) (23) 
for all q sufficiently close to p, q < p. 
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Proof, (i) follows directly from the variational characterization of \ p and rt2J), since 



liv*v,gllp _ _ 

ii ii n ii iiti ^(1 



In order to prove (ii) we note from Theorem [8] that 



lim ||« "'' 



„ i lim ||uJ| p 
y ii q ii y ii p 



since M g converges uniformly to e p when q — > p . Thus, since 



Aq = fi- 



ll) 



J n,q\\ q 



\Ur, 



u 



\m \\ p ~ q \\u \\ P 

Felloe H^yllp 



OWq 



IP' 
yllp 



(24) 



(25) 



we obtain 



lim Aq 

y-s>p- 



lim fi q 

y-s>p- 



lim 

?->p - II w 



9lly 



Ap. 



yllp 



Now we prove error estimate (23). It follows from (i) and (22) that 

< ^ < A ? . 

Hence, 

< max {(Ap - /i 9 ) , (A, - A p )} < A g - /x ff . 
Thus, in order to prove (iii) we need only to bound A q — fi q . It follows from (25) that 

\1 



Ag fiq fiq 



u. 



1\\q 



Vq- 



) dx 



In u qd x 



Therefore, 



A g - fig ^ A ; 



la H ~ 4) dx 



In u Q dx 



A, 



max (t q - t p ) 



\p\Q\ fq\^p-q 



dx 



I n U P q dx \p 

X P \Q\ . 
<^^(p-q). 



V 



In u qdx 

Taking into account (24), there exists R > such that j n u p dx ^ R for all q near to p~. Thus, 

A p |fi| 



II ^yiiy 

< HT, Hp ~Hq^ 



R 



(p-q) = K(p-q). 
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5 Some numerical results 



5.1 Unit Balls 

In this section we present some numerical results in the unit ball of dimensions N = 2,3,4 
applying Algorithm 1, since in this case k p is explicitly known. Computations were performed on 
a Windows 7/ i5 - 4.0 GHz platform, using the GCC compiler. The numerical approximations 
for the first eigenpair were obtained choosing \x = k p and taking q = p — 0.01. Thus, according 



to (13) 



H = k p = N 



p 



p — 1 



We recall from (15) that for the unit ball the functions in the sequence of iterates are radially 
(r = \x\) given by 

1 9 1 

v n +i (r) = J (J {])Y 1 k P v n(s) q - l d.^J " d9, with v {r) = (l - |r|^) . 
Thus, starting with the function 

Vo(r) = ^1 — Irl^ 1 ^ 
we have implemented the sequence of iterates 

v n+l (r) = j^[ [[ 'vM^dsy 1 d9 (26) 



which, after normalized by the sup norm, should be close to the normalized first eigenfunction 
The first eigenvalue was approximated by the sequence 



v 



N ( p 



"Moo XJ 



p-1 



given by ( 19 ). 



In order to compute sequence (26) we mixed the composite Simpson and trapezoidal methods 



on a 101 points mesh for computation of the associated integrals. We adopted 

< 10~ 9 (27) 



|^n+l ^-'n||oo -if-) — !) 



as a stopping criterion. 

At Table [TJ the results for the first eigenvalue of the p-Laplacian for values of p ranging from 
1.1 to 4.0 for the unit balls of dimensions N = 2, 3 and 4 are displayed and truncated at the 
fourth decimal place. The results compare very well with the ones presented in [TD] up to the 
second decimal digit. 
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Table 1: First eigenvalues of the p-Laplacian on the unit ball. 



V 




iV = 2 


N = 3 


N = 4 


V 




N = 2 


N 


= 3 


N 


= 4 


1. 


.1 


2, 


.5666 


3.8665 


5.17607 


2 


.6 


8.08856 


14, 


.9747 


23 


.8345 


1. 


.2 


2. 


,9601 


4.5026 


6.0797 


2 


,7 


8.50354 


15 


.9521 


25 


.672 


1. 


.3 


3. 


.3182 


5.1098 


6.97306 


2 


,8 


8.92654 


16 


.9646 


27 


.6004 


1. 


.4 


3. 


,6637 


5.71889 


7.89478 


2 


,9 


9.35759 


18 


.013 


29 


.6225 


1. 


.5 


4. 


,0053 


6.3419 


8.86046 


3 


.0 


9.79673 


19 


.0977 


31 


.7409 


1. 


.6 


4. 


,3477 


6.98495 


9.87865 


3 


.1 


10.244 


20 


.2194 


33 


.9581 


1. 


.7 


4. 


,6932 


7.65165 


10.955 


3 


.2 


10.6994 


21, 


.3785 


36 


.2769 


1. 


.8 


5. 


,0434 


8.34438 


12.094 


3 


.3 


11.163 


22 


.5755 


38, 


.6999 


1. 


.9 


5. 


,3993 


9.06487 


13.2991 


3 


.4 


11.6347 


23 


.8111 


41, 


.2298 


2. 


.0 


5. 


,7616 


9.81443 


14.5735 


3 


,5 


12.1146 


25 


.0856 


43 


.8694 


2. 


.1 


6. 


,1308 


10.5942 


15.9202 


3 


,6 


12.6027 


26 


.3997 


46 


.6213 


2. 


.2 


6. 


,5071 


11.405 


17.3421 


3 


,7 


13.099 


27 


.7539 


49 


.4884 


2. 


.3 


6. 


,8909 


12.2478 


18.8418 


3 


,8 


13.6034 


29 


.1486 


52 


.4734 


2. 


.4 


7. 


2823 


13.1232 


20.422 


3 


,9 


14.1161 


30 


.5844 


55 


.5792 


2. 


.5 


7. 


,6815 


14.0319 


22.0855 


4 


,0 


14.6369 


32 


.0618 


58, 


.8085 



Graphs of some eigenfunctions generated by the inverse iteration of sublinear supersolutions 
are presented in Figures [IJ [2] and [3] for iV = 2, 3 and 4, respectively. In these graphs it is possible 
to observe the asymptotic behavior of the L°°-normalized eigenfunctions e p with respect to p for 
both cases: p — > l - and p — > oo. The eigenfunctions e p converge to the characteristic function of 
the ball as p — > 1 + (see [35]). On the other hand (see [33]), as p — > oo these functions converge 
to the distance function to the boundary, which in this case is 1 — \x\ . 

Figure [4] illustrates the log concavity of the eigenfunctions e p . Note from Figures [TJ [2] and [3] 
that each eigenfunction e p seems to be convex near the boundary (r = 1). However, log(e p ) is 
surely concave for convex domains, as proved in [43J. 

In Figure 5 we see that ^/Xp approaches 1 as p increases, which is coherent with the following 



known asymptotic behavior (see [33]): lim ^/\, = 1/R where R is the inradius of the domain 



(that is, the radius of the largest ball that lies within the domain). Moreover, one observes from 
Table [TJ that X p approaches the value N of the dimension as p — > l~. It is known (see [35]) that 
X p tends to N / R, if the domain is a ball of radius R. 

Finally, for comparison we show in Figure [6] graphs of p versus X p obtained in two ways: one of 
them through the method proposed in [TO] which is directly based on the inverse power method 
(IPM), while the other is the inverse iteration of sublinear supersolutions (IISS) as developed in 
the present work. 
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Figure 1: Radial profiles of the first eigenfunction for unit ball when N = 2, p = 1.1, 1.2, 1.3, 1.4 
(left) and p = 2.5,3.0,3.5,4.0 (right). 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



Figure 2: Radial profiles of the first eigenfunction for unit ball when N = 3, p = 1.1, 1.2, 1.3, 1.4 
(left) and p = 2.5, 3.0, 3.5, 4.0 (right). 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



Figure 3: Radial profiles of the first eigenfunction for unit ball when N = 4, p = 1.1, 1.2, 1.3, 1.4 
(left) and p = 2.5, 3.0, 3.5, 4.0 (right). 
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Figure 4: Graphs of log(e p ) versus p for the iV-dimensional unit ball, and N = 2 (above), N = 3 
(center), N = 4 (below), and p = 1.5,2.0,2.5,3.0,3.5,4.0. 
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1.18 




"5 11 1 5 2 2 5 i 1 




SO 1 1 5 2 2 5 5 



I .2 2 




5 1 1 5 2 2 5 5 



Figure 5: Graphs of ^/\> versus p for the N- dimensional unit ball, and N = 2 (above), 
(center), N = 4 (below), from p = 50 to p = 290, step 10. 
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Figure 6: Graphs of X p versus p for the N- dimensional unit ball and N = 2 (left), N = 3 (center), 
N = 4 (right). 

5.2 Square, Cube and Torus 

To compute eigenvalues on more general domains, we use a p- version finite element discretization 
on unstructured hexahedral meshes. The discrete equations are solved with PETSc [H] using 
a Newton-Krylov method in which a matrix associated with a lowest-order discretization is 
assembled for preconditioning, while the high-order operator is applied in unassembled form (see 
[T8] for details). For these more complicated domains we apply Algorithm 2. This produces the 
system 

-V- ((e 2 + |V0 m+1 | 2 )^V0 m+1 ) = ( 

where e = 1CT 5 is the regularization used to avoid the singularity or degeneracy at V0 = 0. The 
initial guess for the Newton iteration is taken to be <f) m which leads to very fast convergence in 
the terminal phase. To solve 2D problems with the 3D discretization, homogeneous Neumann 
boundary conditions are imposed on both faces in the z direction. The source code is publicly 
available from https : / / github . com/ j edbrown/ dohp , 

Table [2] shows computed eigenvalues for the unit square and unit cube. These solutions 
were computed using Q5 elements and are as accurate as double precision rounding error for 
the smooth solutions in the p = 2 case. The accuracy of the discretization for a given smooth 
solution has been verified to be essentially independent of p using the method of manufactured 
solutions. This indicates that the primary source of error in Table [2] is interpolation error, as 
usual for finite element methods. 

Figure [7] shows computed eigenfunctions for p = 1.2 and p = 5 on a torus. The unstructured 
hexahedral mesh was created with CUBIT version 13.0 [13] using the commands 

create torus major radius 1 minor radius 0.4 
webcut volume all with plane xplane offset 
mesh volume 1 2 

and a Q2 discretization was used. The computed eigenvalues are \ p = 7.800846 for p = 1.2 and 
A p = 2064.08 for p = 5. 
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Figure 7: Eigenf unctions of the £>-Laplacian for p — 1.2 (top) and p = 5 (bottom) computed on 
a torus with major radius 1 and minor radius 0.4. 
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Table 2: Eigenvalues of the p-Laplacian on the unit square and cube. The 2D results use a 10 x 10 
mesh of Q§ elements, the 3D results use a 6 x 6 x 6 mesh of Q§ elements. For p = 2, the exact 
solutions 27r 2 and 3ir 2 are available so we show the error in the Reference column; these cases are 
denoted by [*]. 

2D 3D 



p 


Computed 


Reference 


Computed 


Reference 


1.2 


6.195550328210643 




8.642315135978254 




1.5 


10.07201415299496 


10.0722 P3] 


14.47791516619582 




1.75 


14.28146165697044 


14.2815 |H] 


20.96672431961172 




2 


19.73920880217817 


5.36xl0" 13 [*] 


29.60881320326431 


3.77xl0- 12 [*] 


2.2 


25.24862830212583 


25.2412 P3J] 


38.51651963302274 




2.5 


35.94868349730170 


35.9493 [H] 


56.19031685699854 




3 


62.75762286200781 


62:7633 [H] 


101.8697977481977 




4 


176.5980821441738 


176.693 Hi] 


306.1647710559179 




5 


463.8206306371868 




849.9777670614186 





Table 3: Convergence of the computed eigenvalue \ m — > X p as q — > p~ for the unit cube using 
a 4 x 4 x 4 mesh with discretization. 

p = 1.5 p = 3 



p- 


• Q 




Xp ^p,q 






Xp ^p,q 




10" 


-l 


13.797661713072 


6.7943x10" 


-l 


96.414190427672 


5.4559 




10" 


-2 


14.405694866696 


7.1393x10- 


-2 


101.31034449471 


5.5975x10- 


-i 


10" 


-3 


14.469912150762 


7.1756x10- 


-3 


101.81397339451 


5.6119x10- 


-2 


10" 


-4 


14.476369813850 


7.1792x10- 


-4 


101.86447907670 


5.6133x10- 


-3 


10" 


-5 


14.477015941408 


7.1796x10- 


-5 


101.86953107554 


5.6135x10- 


-4 


io- 


-6 


14.477080557778 


7.1796x10- 


-6 


101.87003628973 


5.6135x10- 


-5 







14.477087737416 






101.87009242480 







Experimental evidence suggests that Algorithm 2 converges with q = p, but we have only 
been able to prove convergence for q < p. It is unknown whether the iteration will break down 
for some domain when q = p, but one can always compute with q < p in which case Theorem 8 
guarantees convergence with an error less than K(p — q) for some positive constant K depending 
only on the domain and p. Table [3] shows numerical evidence of this result and quantifies K for 
the unit cube with p = 1.5 and p = 3. 

In practice, the total computational cost to solve the eigenvalue problem is about twice that 
of only solving the torsion creep problem. Table [4] shows the convergence of inverse iteration 
when the Newton iteration at each step is started using the solution at the last iteration. The 
initial guess for the torsion creep problem is zero, which leads to a difficult nonlinear solve. 
The Newton iteration is guarded by a cubic backtracking line search which is sufficient in this 
case; a parameter continuation or grid sequencing is more robust. The torsion creep problem is 
significantly easier to solve for less extreme values of p or for larger values of the regularization e. 
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Table 4: Convergence rate for inverse iteration applied to the torus with p — 1.2 and p — 5. Each 
nonlinear solve is converged to a relative tolerance of 1CT 8 . 

p — 1.2 p = 5 



Newton its. 




Newton its. 




37 


torsion 


18 


torsion 


5 


7.7670871 


6 


1628.81 


4 


7.7965212 


4 


1975.40 


3 


7.8003037 


4 


2043.11 


3 


7.8007802 


3 


2057.15 


2 


7.8008389 


3 


2061.17 


2 


7.8008456 


3 


2062.74 


2 


7.8008462 


3 


2063.35 






3 


2063.38 






3 


2063.98 






3 


2064.08 



After the torsion creep problem has been solved, a line search is no longer necessary and accurate 
estimates of the eigenvalue can be obtained in a few more Newton iterations. 
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